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Abstract 

We present a preferential attachment growth model to obtain the distribution P{K) of number 
of units K in the classes which may represent business firms or other socio-economic entities. We 
found that P(K) is described in its central part by a power law with an exponent ip = 2 + b/{\ — b) 
which depends on the probability of entry of new classes, b. In a particular problem of city 
population this distribution is equivalent to the well known Zipf law. In the absence of the new 
classes entry, the distribution P(K) is exponential. Using analytical form of P(K) and assuming 
proportional growth for units, we derive P(g), the distribution of business firm growth rates. The 
model predicts that P(g) has a Laplacian cusp in the central part and asymptotic power-law tails 
with an exponent £ = 3. We test the analytical expressions derived using heuristic arguments by 
simulations. The model might also explain the size-variance relationship of the firm growth rates. 
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I. INTRODUCTION 

Gibrat , building upon the work of the astronomer Kapteyn |3j , assumed the expected 
value of the growth rate of a business firm's size to be proportional to the current size of the 
firm, which is called "Law of Proportionate Effect" |5j. Several models of proportional 
growth have been subsequently introduced in economics in order to explain the growth of 
business firms jgl Q, 0| • Simon and co-authors extended Gibrat's model by 

introducing an entry process according to which the number of firms rise over time. In 
Simon's framework, the market consists of a sequence of many independent "opportunities" 
which arise over time, each of size unity. Models in this tradition have been challenged by 



many researchers 
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12 1 who found that the firm growth distribution is not 



Gaussian but displays a tent shape. 

Here we introduce a mathematical framework that provides an unifying explanation for 
the growth of business firms based on the number and size distribution of their elementary 
constituent components HyHHHHHQ. Specify we present a m odel of 
proportional growth in both the number of units and their size and we draw some general 



implications on the mechanisms which sustain business firm growth 
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According to the model, the probability density function (PDF) of growth rates, P(g) is 

13 n 



Laplace [29] in the center [l3j] with power law tails 3(J,|31| decaying as where ( = 3. 

Two key sets of assumptions in the model are described in subsections A (the number 
of units K in a class grows in proportion to the existing number of units) and B (the 
size of each unit fluctuates in proportion to its size). Our goal is to first find P(K), the 
probability distribution of the number of units in the classes at large t, and then find 
P(g)usingtheconvolutionofP(K) and the conditional distribution of the class growth rates 
P(g\K), which for large K converges to a Gaussian.. 



II. ANALYTICAL RESULTS 

A. The Proportional Growth of Number of Units 

The first set of assumptions is: 



(Al) Each class a consists of K a (t) number of units. At time t = 0, there are N(0) classes 
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consisting of n(0) total number of units. The initial average number of units in a class 
is thus n(0)/JV(0). 

(A2) At each time step a new unit is created. Thus the number of units at time t is 
n(t) = n(0)+t. 

(A3) With birth probability 6, this new unit is assigned to a new class, so that the average 
number of classes at time t is N(t) = N(0) + bt. 

(A4) With probability 1 — b, a new unit is assigned to an existing class a with probability 
P a = (1 - b)K a (t)/n(t), so K a (t + 1) = K a (t) + 1. 

This model can be generalized to the case when the units are born at any unit of time t' 
with probability /i, die with probability A, and in addition a new class consisting of one unit 
can be created with probability b' [32]. This model can be reduced to the present model if 
one introduce time t = t'(fi — A + b') and probability b = b' /(// — A + b'). 

Our goal is to find P(K), the probability distribution of the number of units in the classes 
at large t. This model in two limiting cases (i) 6 = 0, K a = 1 (a = 1, 2 . . . N(0)) and (ii) 



01= 1, ra(0) = 1 has exact analytical solutions P(K) = N(0)/t(t/(t + N(0))) K (l + 



3j, El and lim P(K) = (1 + b)T(K)T(2 + b)/T{K + 2 + 6) |35( respectively, In 



6 ^ 0, N 

o{i/t)) 

t— >oo 

general, an exact analytical solution of this problem cannot be presented in a simple close 
form. Accordingly, we seek for an approximate mean-field type (3^ solution which can be 
expressed in simple integrals and even in elementary functions in some limiting cases. First 
we will present a known solution of the preferential attachment model in the absence of the 
influx of new classes (37 1: 

PoUK) = ^x^l « eM-K/K(tW + 0(t~% (1) 

where A = 1 — l/K(t) and K(t) = [n(0) + t]/N(0) is the average number of units in the old 
classes at time t. Note that the form of the distribution of units in the old classes remains 
unchanged even in the presence of the new classes, whose creation does not change the 
preferential attachment mechanism of the old classes and affects only the functional form of 
K(t). 

Now we will treat the problem in the presence of the influx of the new classes. Assume 
that at the beginning there are N(0) classes with n(0) units. Because at every time step, 
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one unit is added to the system and a new class is added with probability b, at moment t 
there are 

n{t)=n{0)+t (2) 

units and approximately 

N(t) = N(0) + bt (3) 

classes, among which there are approximately bt new classes with n new units and iV(0) old 
classes with n i d units, such that 

n Q id + n new = n(0) + t. (4) 

Because of the preferential attachment assumption (A4), we can write, neglecting fluctua- 
tions j3(| and assuming that t, n i d , and n new are continuous variables: 

d<Tl new . « f-neu) / r \ 
: = + (l — OJ— 7—, , (5) 

dt v J n(0) +t K J 

dn id . . n id , > 

~w = {1 - b) ^on-f (6) 

Solving the second differential equation and taking into account initial condition n o ^(0) = 
n(0), we obtain n id(t) = (n(0) + n(0) b . Analogously, the number of units at time t in 
the classes existing at time to is 

n e (t ,t) = (n(0)+t) l ~ b (n(0)+t ) b (7) 

where the subscript 'e' means "existing". Accordingly, the average number of units in old 
classes is 

K(+ , n old {t) (n(0)+ty- b b 
Thus according to Eq. ((TJ), the distribution of units in the old classes is 



p r A^(0) / KN(0) x 

old{ ' ^ (n(0)+ty- b n(0) b 6XP \(n(0)+ty- b n{0) b >' 



and the contribution of the old classes to the distribution of all classes is 

P old {K) = P old (K)N(0)/(N(0) + bt). (10) 

The number of units in the classes that appear at to is b dt and the number of these classes is 
b dt. Because the probability that a class captures a new unit is proportional to the number 



of units it has already gotten at time t, the number of units in the classes that appear at 
time to is 

n new {t , t) = n e (t , t)bdt/[n(0) + t }. (11) 
The average number of units in these classes is 

K(t ,t) = n new (t ,t)/bdt = (n(0)+t) 1 - b /(n(0)+t ) 1 - b . (12) 

Assuming that the distribution of units in these classes is given by a continuous approxima- 
tion we have 

P new {K, t ) ~ "T777 it exp {-K/K{t , £)) . (13) 
K(t ,t) 

Thus, their contribution to the total distribution is 

b dtn 1 . i j . . 

The contribution of all new classes to the distribution P{K) is 

P ^ K) ™ nwth t mj) exp ( - K/K{t °' 4)) dt °- (14 » 

If we let y = K/K(t ,t) then P new (K) = P new (K)bt/(N(0) + bt) where 

P„UK) « ^)A±i K f „ A (15) 

and the low limit of integration, K' is given by 

K , = K («2L_Y- b (16) 



v n(0) + t 

Finally the distribution of units in all classes is given by 

Now we investigate the asymptotic behavior of the distribution in Eq. ()15j) and show that 
it can be described by the Pareto power law tail with an exponential cut-off. 
1. At fixed K when t — > oo, we have K' — ► 0, thus 

i 







1-6 



oo 



K 



K~ L -—K (18) 



As if — > oo, P new (K) converges to a finite value: 

P new (K) = K-'-\ T (l + J-^) . (19) 

Thus for large if ^> 1, but such that if' <C 1 or if <C (l+i/n(0) 1-6 , we have an approximate 
power-law behavior: 

P new (K) ~ if ^, (20) 

where y2 = 2 + 6/(1 - 6) > 2. 
As if ^ 0, 

1 f 1 ^ if( 1+ T^) 1 
iW*0 = — ~T if™^ t— = — -T- (21) 

v ; 1-6 i + _i_2-6 v; 

2. At fixed t when if — > oo, we use the partial integration to evaluate the incomplete T 
function: 

/■OO /"OO 

-v y a dy = -e- y y a \™ + a e~ y y a ~ x dy « e - * x a . 

«/ a; 

Therefore, from Eq. (|15j) we obtain 

P_(if) - A if^- 1 r e- ^ dy, 

n(o) + u i-b J K (^r b 

n{0) ' I -p (-K , (22) 



iV(0) + 6* 1 - b if \ \n(0)+t 

which always decays faster than Eq. (JUJ) because n(0) > N(0) and there is an additional 
factor if -1 in front of the exponential. Thus the behavior of the distribution of all classes 
is dominated for large if by the exponential decay of the distribution of units in the old 
classes. 

Note that Eq. Q and Eq. ()15j) are not exact solutions but continuous approximations 
which assume if is a real number. This approximation produces the most serious discrepancy 
for small if. To test this approximation, we perform numerical simulations of the model for 
b = 0.1, N(0) = n(0) = 10000 and t = 400000. The results are presented in Figffl While 
the agreement is excellent for large if, Eq. (|15|) significantly underestimates the value of 
P new (K) for if = 1 and if = 2. Note that in reality the power- law behavior of P new (K) 
extends into the region of very small if. 
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B. The Proportional Growth of Size of Units 



The second set of assumptions of the model is: 

(Bl) At time t, each class a has K a (t) units of size i = 1, 2, ...K a (t) where K a and 

£i > are independent random variables taken from the distributions P(K a ) and 
Pg(£i) respectively. P(K a ) is defined by Eq. (fT7|) and P^(Ci) is a given distribution 
with finite mean and standard deviation and ln£j has finite mean /i£ = (ln£j) and 
variance = ((ln£j) 2 ) — The size of a class is defined as S a (t) = 

(B2) At time t + 1, the size of each unit is decreased or increased by a random factor 
rjiit) > so that 

+ = (23) 

where ^j(t), the growth rate of unit 2, is independent random variable taken from a 
distribution P v (r]i), which has finite mean and standard deviation. We also assume 
that In 7^ has finite mean /i v = (lnr^) and variance = ((ln^) 2 ) — ft 2 . 

Let us assume that due to the Gibrat process, both the size and growth of units (£j and 
r/i respectively) are distributed lognormally 

p(6) = -jL= i exp (-(ln& - m ? ) 2 /2^) , (24) 
I i 

p(Vi) = 7== — exp (-Qnrti - m ri ) 2 /2V v ) . (25) 

If units grow according to a multiplicative process, the size of units £■ = ^rji is distributed 
lognormally with = + V v and = + m v . 

The n th moment of the variable x distributed lognormally is given by 

f°° 1 x n 

^x(n) = / —= — dx exp (— (lnx — m) 2 /2V) = exp (nm x + n 2 V x /2) . (26) 
Jo \'2ttV x 

Thus, its mean is \i x = fi x (l) = exp(m x + V x /2) and its variance is a 2 = /i^) — /^(l) 2 = 
/ i ;c (l) 2 (exp(V;)-l). 

Let us now find the distribution of g growth rate of classes. It is defined as 

^ln^±ll = lnX:C'-lnX:^ (27) 
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Here we neglect the influx of new units, so K a = K a (t + 1) = K a (t). 

The resulting distribution of thegrowth rates of all classes is determined by 

oo 

P(g)^J2 P ( K ) P (9\n (28) 

K=l 

where P{K) is the distribution of the number of units in the classes, computed in the 
previous stage of the model and P(g\K) is the conditional distribution of growth rates of 
classes with given number of units determined by the distribution P^(0 and -P^(^). 

Now our goal is to find an analytical approximation for P(g\K). According to the central 
limit theorem, the sum of K independent random variables with mean = fi^(l) and finite 
variance <r| is 

K 

J2^i = K^ + VKu K , (29) 
i=i 

where vk is the random variable with the distribution converging to Gaussian 

lim P(u K ) - — L= exp (-4/2^|) • (30) 

Accordingly, we can replace ln(5^ i=1 by its Tailor's expansion In K + In fi^ + vk / (^VK) , 
neglecting the terms of order K~ l . Because In \i^ = + V v /2 and In [i? = ln/ig + In fi v we 
have 



g = \n S(t + 1) - In S(t) = \n(K/i e ) + H K ^) ~ 



= m ^ + '_mp^. (31) 

For large K the last term in Eq. (|31jl is the difference of two Gaussian variables and that 
is a Gaussian variable itself. Thus for large K, g converges to a Gaussian with mean, 
m = + V v /2, and certain standard deviation which we must find. 
In order to do this, we rewrite 



y/Kfi? K 11? 

and 
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Thus 

= m v + ^+ ^^-^ . (32) 
2 K n^i 

Since fi^ = /x^/i,,, the average of each term in the sum is fi^i — /i£ fi v = 0. The variance of 
each term in the sum is rji) 2 ) — (2£ 2 rji fi v ) + (£ 2 /£) where ^rji, Q^r\i and £ 2 are all lognormal 
independent random variables. Particularly, (£j?7i) 2 is lognormal with V = 4V V + 4Vg and 
m = 2m v + 2m^; ^??7i is lognormal with V = 4V^ + and m = 2m^ + m v ; £ 2 is lognormal 
with V = 4Vg and m = 2mg. Using Eq. (J2f 



((tiVi) 2 ) = exp(2m„ + 2m ? + 2V V + 2V i ), (33a) 
(gin) = exp(m r , + 2m € + 2V^ + Kj/2), (33b) 
(C 2 )=exp(2m ? + 2^). (33c) 

Collecting all terms in Eqs. (|33aH33c|) together and using Eq. (J3*2j) we can find the variance 
of 5-: 

2 K exp(2m c + 2V C + 2m v + V v )(exp(V v ) - 1) 
a ~ K 2 exp(2m € + V^ + 2m v + V v ) ' 

= l e xp(^)(exp(K ? )-l). (34) 
Therefore, for large K, g has a Gaussian distribution 

where 

m = m v + V v /2 (36) 

and 

V = Ka 2 = exp(Vz)(exp(V v ) - 1). (37) 

Note, that the convergence of the sum of lognormals to the Gaussian given by Eq. (J29|) 
is a very slow process, achieving reasonable accuracy only for K /^(2) ~ exp(2V^). For 
a pharmaceutical database [3j|, we have = 5.13, = 3.44, V v = 0.36, and m v = 0.16. 
Accordingly, we can expect convergence only when K ^> 3 • 10 4 . Figur^21 demonstrates the 
convergence of the normalized variance Ka 2 (K) and mean m(K) of g to the theoretical 



limits given by Eqs. (jHBj) and (J37|) respectively: V = 73.24 and m = 0.196. In both 
cases, the discrepancy between the limiting values and the actual values decreases as l/\/K. 
Interestingly, Eq. (f33|) predicts cr(K) ~ K" 13 , where 3 = 1/2. This value is much larger 
than the empirical value 3 ~ 0.2 observed for the size-variance relationships of various 
socio-economic entities [3, 13, Q ■ However, the slow convergence of V(K)K suggests 
that for quite a wide range of K < 1000, a(K) ~ K~ 2 and only at K > 10 4 there is a 
crossover to the theoretical value 3 = 0.5, (Fig. |3J). Finally, the simulated distribution of 
P(g\K) has tent-shape wings which develop as K increases (Fig. HJ). This feature of the 
model growth rates may explain the abundance of the tent-shaped wings of the growth rates 
of various systems in nature. The most drastic discrepancy between the Gaussian shape 
and the simulated distribution P(g\K) can be seen when K m 1000 and than it starts to 
decrease slowly, and remains visible even for K = 10 6 . 

Nevertheless, in order to obtain close form approximations for the growth rate, we will 
use the Gaussian approximation (}35|) for P(g\K). The distribution of the growth rate of the 
old classes can be found by Eq. (|2*K|) . In order to find a close form approximation, we replace 
the summation in Eq. (|2~£j) by integration and replace the distributions P(K) by Eq. (jHJ) and 
P(g\K) by the Eq. (J33j) . Assuming m = 0, we have 

1 f°° 1 K g 2 K 



= } mr 1 (38) 

2V2V \ 2V y J ' 1 ; 

where K(t) is the average number of units in the old classes (see Eq. (jSJ)). This distribution 
decays as 1/g 3 and thus does not have a finite variance. In spite of drastic assumptions that 
we make, Eq. (J3~%j) correctly predicts the shape of the convolution P id(g). Figure El shows 
the comparison of the simulation of the growth rates in the system with the exponential 
distribution of units P(K) with K(t) = 2 15 and the same empirical parameters of the 
unit size and growth distributions as before. The parameter of the analytical distribution 
characterizing its width (variance does not exist), must be taken V = 33 which is much 
smaller than the analytical prediction V = 73.23. This is not surprising, since for K = 2 15 
Ka 2 (K) = 50 (see Fig. Eb). Moreover, since we are dealing with the average a 2 (K)K for 
K < 2 15 , we can expect V < 50. Nevertheless the nature of the power-law wings decaying 
as 1/g 3 is reproduced very well. 



10 



For the new classes, when t — > oo the distribution of number of units is approximated by 



1 '* 



p ncw (K) « Y^i K ^~ b I y~ be ~ Vd y- ( 39 ) 







Again replacing summation in Eq. (|28|) in the text by integration and P(g\K) by Eq. (|35|) 
and after the switching the order of integration we have: 



J | poo poo 

Pne W (g)~Y^l^fJ eM^y^dy J ex V {-g 2 K/2V)K { -^-^ dK. 

As g — > oo, we can evaluate the second integral in Eq. (jlpj) by partial integration: 

1 f°° 1 2V ill 
Pnew(g) ~ -, — r / j== — y~^~* y^~ b exp(-y) exp(-y^ 2 / 2 ^) 
1 - o io V2irV g 2 



(40) 



1 1 2V 1 1 

X 7T ~ 



(41) 



1 - b ^/2W 9 2 ^g 2 /2V + 1 £ 3 ' 

We can compute the first derivative of the distribution (|4T)|) by differentiating the inte- 
grand in the second integral with respect to g. The second integral converges as y — > 0, 
and we find the behavior of the derivative for g —>■ by the substitution x = Kg 2 /(2V). 
As g — > 0, the derivative behaves as g g 2 \~(z/2)+i/{'i--b)\ ~ g 2b K l - b )^ which means that the 
function itself behaves as C2 — Ci|<5>| 2b /( 1_b ) +1 , where C2 and G\ are positive constants. For 
small b this behavior is similar to the behavior of a Laplace distribution with variance V: 
exp(-V2\g\/VV)/V2V = 1/V2V -\g\/V. 

When b — > 0, Eq. ()40|) can be expressed in elementary functions: 

1 POO pK 

PneM\b->o « -f= \ K-^ 2 eM~Kg 2 /2V)dK / exp{-y)ydy, 
V2ttV Jo Jo 

l/l 2 

~ Tff \ v/l + gi/2 V + \g\/V2V + vV/2 V+T 

Simplifying we find the main result: 

2V 

PneMlb^O ~ , ~- — , ~- (42) 

which behaves for g — > as \ j\/2V — \g\jV and for g — > 00 as V/ (2g 3 ). Thus the distribution 
is well approximated by a Laplace distribution in the body with power-law tails. Because of 
the discrete nature of the distribution of the number of units, when g ^> \J2V the behavior 
for g — > 00 is dominated by exp(— g 2 /2V). 
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In Fig. we compare the distributions given by Eq. (|38jh the mean field approximation 
Eq. ()4Uj) for b = 0.1 and Eq. (}4"2*|) for b 0. We find that all three distributions have very 
similar tent shape behavior in the central part. In Fig. |Hb we also compare the distribution 
Eq. (J42j) with its asymptotic behaviors for g — > (Laplace cusp) and g —>■ oo (power law), 
and find the crossover region between these two regimes. 

III. CONCLUSIONS 

The analytical solution of this model can be obtained only for certain limiting cases but 
a numerical solution can be easily computed for any set of assumptions. We investigate the 
model numerically and analytically (see and find: 

(1) In the presence of the influx of new classes (b > 0), the distribution of units converges 
for t -> oo to a power law P{K) ~ K^, ip = 2 + 6/(1 - b) > 2. Note that this 
behavior of the power-law probability density function leads to the power law rank- 
order distribution where rank of a class R is related to the number of its units K as 



Thus K ~ i where ( = l/(<p — l) = l — b< 1, which leads in the limit b — > to the 
celebrated Zipf's lawj^J for cities populations, K ~ 1/R. Note that this equation can 
be derived for our model using elementary considerations. Indeed, due to proportional 
growth the rank of a class, R, is proportional to the time of its creation t . The 
number of units n(to) existing at time to is also proportional to to and thus also 
proportional to R. According to the proportional growth, the ratio of the number of 
units in this class to the number of units in the classes existed at time to is constant: 
K(to,t)/n e (to,t) = l/n(to). If we assume that the amount of units in the classes, 
created after t can be neglected since the influx of new classes b is small, we can 
approximate n e (to,t) « n(t) ~ t. Thus for large t, n e (t ,t) is independent of to and 
hence K(to,t) ~ 1/R. If we do not neglect the influx of new classes, Eq. (JJJ) gives 
n e (t ,t) ~ 4, hence K{t ,t) ~ l/R l - b . 

(2) The conditional distribution of the logarithmic growth rates P(g\K) for the classes 
consisting of a fixed number K of units converges to a Gaussian distribution 




(43) 
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for K — y oo. Thus the width of this distribution, <r(K), decreases as l/K 13 , with 
(3 = 1/2. Note that due to slow convergence of lognormals to the Gaussian in case 
of wide lognormal distribution of unit sizes = 5.13, computed from the empirical 
data j^J, we have (3 = 0.2 for relatively small classes. This result is consistent with 
the observation that large firms with many production units fluctuate less than small 



firms 



7 



18 



22 



381 ] . Interestingly, in case of large V£, P(g\K) converges to the Gaus- 
sian in the central interval which grows with K, but outside this interval it develops 
tent-shape wings, which are becoming increasingly wider, as K — > oo. However, they 
remain limited by the distibution of the logarithmic growth rates of the units, P^ilnrj). 

(3) For g 3> V v , the distribution P(g) coincides with the distribution of the logarithms of 
the growth rates of the units: 

P(g)^P v (lnr]). (44) 

In the case of power law distribution P(K) ~ K"^ which dramatically increases for 
K — > 1, the distribution P(g) is dominated by the growth rates of classes consisting 
of a single unit K — 1, thus the distribution P (q) practically coincides with P rj (hirj) 
for all g. Indeed, empirical observations of Ref. |3j| confirm this result. 

(4) If the distribution P{K) ~ K~ v , ip > 2 for K — > oo, as happens in the presence of the 
influx of new units 67^0, P(g) = C\ — C2\g\ 2ip ~ 3 , for g — > which in the limiting case 
b — ► 0, ip —> 2 gives the cusp P(g) ~ C\ — C^g\ [C\ and C2 are positive constants), 
similar to the behavior of the Laplace distribution Ph{g) ~ exp(— |p|C 2 ) for g — > 0. 

(5) If the distribution P{K) weakly depends on K for K — > 1, the distribution of P(g) 



can be approximated by a power law of g: P(g) ~ (? _3 in wide range y/V g /K(t) <S 
5» \A^, where is the average number of units in a class. This case is realized 
for b = 0, t — > 00 when the distribution of P(K) is dominated by the exponential 
distribution and i^(t) — >• 00 as defined by Eq. ([!}. In this particular case, P(g) for 
g can be approximated by Eq. (JHEJ) 

(6) In the case in which the distribution P{K) is not dominated by one-unit classes but 
for K —>■ 00 behaves as a power law, which is the result of the mean field solution for 
our model when t — > 00, the resulting distribution -P(g) has three regimes, P(g) ~ 
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C\ — C2\g\ 2ip ~ z for small g, P(g) ~ g~ 3 for intermediate g, and P{g) ~ P(\nr]) for 
5> — > oo. The approximate solution of -P(g) in this case is given by Eq. (j4*Uj) For 
b 7^ Eq. (J4(J|) can not be expressed in elementary functions. In the b — > case, 
Eq. (J4(J|) yields the main result Eq. (|42|) . which combines the Laplace cusp for g — > 
and the power law decay for g — > oo. Note that due to replacement of summation by 
integration in Eq. (|28jl . the approximation Eq. (|42|) holds only for g < \fV~ r] . 

In conclusion we want to emphasize that although the derivations of the distributions (|38jh 
(|40|) . and (}4*2"j) are not rigorous they satisfactory reproduce the shape of empirical data, 
especially the 1/g 3 behavior of the wings of the distribution of the growth rates and the 
sharp cusp near the center. 
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FIG. 1: Comparison of the distributions P(K) for the new and old classes obtained by numerical 
simulations of the model with the predictions of Eq. ()14|) and Eq. ()10|) respectively. For large 
K the agreement is excellent. The discrepancy exists only for P new at small K, e.g. Eq. 1)141) 
significantly underestimates the P new (l) and P new (2). 
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FIG. 2: Convergence of the parameters of the simulated P(g\K) to the values, which follow from 
the central limit theorem: (a) the mean m(k) and (b) the normalized variance Ka 2 (K). In both 
cases the speed of convergence is 1/VK as can be seen from the straight line fits versus \j\f~K 
with the intercepts equal to the analytical values m = 0.196 and V = 73.24, respectively. The 
parameters of the simulations = 5.13 = 3.44, Vy = 0.36, and = 0.016 are taken from the 
empirical analysis of the pharmaceutical data base 39] . 
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FIG. 3: Crossover of the size-variance relationship <r(K) from K 2 for small K to IT ' 5 for large 
K. The parameters of the simulations are the same as in Fig[2J 

10° 
10"' 

ST 

ra 10^ 

10"' 
10"* 

-10 -5 5 10 

(g-m)/o 

FIG. 4: Convergence of the shape of the distribution of P(g\K) found in simulations to limiting 
Gaussian. One can see the developments of the tent-shape wings as K grows. The parameters of 
the simulations are the same as in Fig|21 
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FIG. 5: Comparison of the shape of the distribution of P(g) for the exponential distribution of 
P(K) = l/{K)exp(-K/{K)) with (K) = 2 15 with the prediction of En. ljSSj) . The parameters of 
the simulation are the same as in Fig. |5J The fitting parameter V = 33 in Eq. ()38|) gives the best 
agreement with the simulation results. One can see a very good convergence to the inverse cubic 
law for the wings. 
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FIG. 6: (a) Comparison of three different approximations for the growth rate PDF, P(g), given by 
Eq. Q38|) . mean field approximation Eq. (|4U|) for b = 0.1 and Eq. (|42j) . Each P(g) shows similar tent 
shape behavior in the central part. We see there is little difference between the three cases, 6 = 
(no entry), b = 0.1 (with entry) and the mean field approximation. This means that entry of new 
classes (6 > 0) does not perceptibly change the shape of P(g). Note that we use K(t)/V g = 2.16 
for Eq. (j38j) and V g = 1 for Eq. (|42|) , (b) The crossover of P(g) given by Eq. (|42|) between the 
Laplace distribution in the center and power law in the tails. For small g, P(g) follows a Laplace 
distribution P(g) ~ exp(— 1<7|), and for large g, P(g) asymptotically follows an inverse cubic power 
law P(g) ~ g- 3 . 
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